#ifndef INT_H
#define INT_H

#include <functional>
#include <vector>
#include <queue>
#include <map>
#include <cmath>
#include <algorithm>
#include <assert.h>
#include <limits>
#include <ostream>

#include "MoreFunctional.h"
#include "SeriesLevin.h"
#include "FFT.h"
#include "IntGKPData10.h"

namespace integ {

#ifdef NDEBUG
#define integ_dbg1 if (false) (*dbgout)
#define integ_dbg2 if (false) (*reg.dbgout)
#define integ_dbg3 if (false) (*tempreg.dbgout)
#else
#define integ_dbg1 if (dbgout) (*dbgout)
#define integ_dbg2 if (reg.dbgout) (*reg.dbgout)
#define integ_dbg3 if (tempreg.dbgout) (*tempreg.dbgout)
#endif

  template <class T> T norm(const T& x) { return x*x; }
  using std::norm;
  template <class T> T real(const T& x) { return x; }
  using std::real;

  //#define COUNTFEVAL  // If defined, then count the number of function evaluations
#ifdef COUNTFEVAL
  int nfeval = 0;
#endif

  enum INTALG {GKP, GKPOSC, OSCINFTWOPI, OSCINFONE};
#define DEFALG GKP
  // GKP = adaptive Gauss-Kronrod-Patterson with bisection for each split
  // GKPOSC = adaptive GKP, with splits at zeros of function
  // OSCINF* = integration of oscillatory functions to infinity
  //       TWOPI = period of oscillation is 2 PI
  //       ONE   = period of oscillation is 1

  template <class T>
    struct IntRegion {
      public:
	IntRegion(const T _a, const T _b, INTALG _alg=DEFALG, std::ostream* _dbgout=0) :
	  a(_a), b(_b), error(0.), area(0), alg(_alg), dbgout(_dbgout) {}
	bool operator<(const IntRegion<T>& r2) const {return error < r2.error;}
	bool operator>(const IntRegion<T>& r2) const {return error > r2.error;}
	void SubDivide(std::vector<IntRegion<T> >* children) {
	  assert(children->size() == 0);
	  if (splitpoints.size() == 0) Bisect();
	  if (splitpoints.size() > 1) std::sort(splitpoints.begin(),splitpoints.end());
	  /*
	     if (a>splitpoints[0] || b < splitpoints.back()) {
	     std::cerr<<"a,b = "<<a<<','<<b<<std::endl;
	     std::cerr<<"splitpoints = ";
	     for(size_t i=0;i<splitpoints.size();i++) { std::cerr<<splitpoints[i]<<"  "; }
	     std::cerr<<std::endl;
	     }
	     */
	  assert(splitpoints[0] >= a);
	  assert(splitpoints.back() <= b);
	  children->push_back(IntRegion<T>(a,splitpoints[0],alg,dbgout));
	  for(uint i=1;i<splitpoints.size();i++) {
	    children->push_back(IntRegion<T>(splitpoints[i-1],splitpoints[i],alg,dbgout));
	  }
	  children->push_back(IntRegion<T>(splitpoints.back(),b,alg,dbgout));
	}
	void Bisect() {splitpoints.push_back((a+b)/2.);}
	void AddSplit(const T x) {splitpoints.push_back(x);}
	uint NSplit() const {return splitpoints.size();}

	void FindOscSplits(const std::map<T,T>& fx)
	{
	  T lastzero = std::min(Left(),Right());
	  T minwidth = fabs(Right()-Left())/20.;
	  typedef typename std::map<T,T>::const_iterator mapit;
	  mapit prevpair = fx.begin();
	  mapit thispair = prevpair; thispair++;
	  while (thispair != fx.end()) {
	    const T& f1 = real(prevpair->second);
	    const T& f2 = real(thispair->second);
	    const T& x1 = prevpair->first;
	    const T& x2 = thispair->first;
	    if (f1 < 0. && f2 >= 0. && x2 > lastzero + minwidth) {
	      T zero = (f2*x1 - f1*x2)/(f2-f1);
	      if (zero - lastzero > minwidth) {
		AddSplit(zero);
		lastzero = zero;
	      }
	    }
	    prevpair = thispair;
	    thispair++;
	  }
	  if (splitpoints.size() == 0) {
	    integ_dbg1<<"No zeros found... bisect\n";
	    Bisect();
	  }
	  else {integ_dbg1<<"Found "<<splitpoints.size()<<" zeros\n";}
	}

	const T& Left() const {return a;}
	const T& Right() const {return b;}
	const T& Err() const {return error;}
	const T& Area() const {return area;}
	INTALG Alg() const {return alg;}
	void SetArea(const T& a, const T& e) {area = a; error = e;}
	void SetAlg(INTALG _alg) {alg = _alg;}

      private:
	T a,b,error,area;
	std::vector<T> splitpoints;
	INTALG alg;

      public:
	std::ostream* dbgout;
    };

  template <class T> T Epsilon() { return std::numeric_limits<T>::epsilon(); }
  template <class T> T MinRep() { return std::numeric_limits<T>::min(); }

/*
  template <> double Epsilon<double>()
  { return 2.2204460492503131e-16; }
  template <> double MinRep<double>()
  { return 2.2250738585072014e-308; }
*/

#ifdef EXTRA_PREC_H
  template <> Quad Epsilon<Quad>()
  { return 3.08148791094e-33; }
  template <> Quad MinRep<Quad>()
  { return 2.2250738585072014e-308; }
#endif

  template <class T>
  inline T rescale_error (T err, const T& resabs, const T& resasc)
  {
    if (resasc != 0. && err != 0.) {
      const T scale = (200. * err / resasc);
      if (scale < 1.) err = resasc * scale * sqrt(scale) ;
      else err = resasc ;
    }
    if (resabs > MinRep<T>() / (50. * Epsilon<T>())) {
      const T min_err = 50. * Epsilon<T>() * resabs;
      if (min_err > err) err = min_err;
    }
    return err ;
  }

  template <class UF>
    int IntGKPNA(const UF& func, IntRegion<typename UF::result_type>& reg,
	const typename UF::result_type epsabs, 
	const typename UF::result_type epsrel,
	std::map<typename UF::result_type,typename UF::result_type>* retfxmap=0)
    // A non-adaptive integration of the function f over the region reg.
    // The algorithm computes first a Gaussian quadrature value
    // then successive Kronrod/Patterson extensions to this result.
    // The functions terminates when the difference between successive
    // approximations (rescaled according to rescale_error) is less than 
    // either epsabs or epsrel * I, where I is the latest estimate of the 
    // integral.
    // The order of the Gauss/Kronron/Patterson scheme is determined
    // by which file is included above.  Currently schemes starting 
    // with order 1 and order 10 are calculated.  There seems to be 
    // little practical difference in the integration times using 
    // the two schemes, so I haven't bothered to calculate any more.
    {
      if (retfxmap) retfxmap->clear();
      typedef typename UF::result_type T;
      const T a = reg.Left();
      const T b = reg.Right();

      const T half_length =  0.5 * (b - a);
      const T abs_half_length = fabs (half_length);
      const T center = 0.5 * (b + a);
      const T f_center = func(center);
#ifdef COUNTFEVAL
      nfeval++;
#endif

      assert(gkp_wb<T>(0).size() == gkp_x<T>(0).size()+1);
      T area1 = gkp_wb<T>(0).back() * f_center;
      std::vector<T> fv1, fv2;
      fv1.reserve(2*gkp_x<T>(0).size()+1);
      fv2.reserve(2*gkp_x<T>(0).size()+1);
      for (uint k=0; k<gkp_x<T>(0).size(); k++) {
	const T abscissa = half_length * gkp_x<T>(0)[k];
	const T fval1 = func(center - abscissa);
	const T fval2 = func(center + abscissa);
	area1 += gkp_wb<T>(0)[k] * (fval1+fval2);
	fv1.push_back(fval1);
	fv2.push_back(fval2);
	if (retfxmap) {
	  (*retfxmap)[center-abscissa] = fval1;
	  (*retfxmap)[center+abscissa] = fval2;
	}
      }
#ifdef COUNTFEVAL
      nfeval+=gkp_x<T>(0).size()*2;
#endif

      integ_dbg2<<"level 0 rule: area = "<<area1<<std::endl;

      T err=0; 
      bool calcabsasc = true;
      T resabs=0., resasc=0.;
      for (int level=1; level<NGKPLEVELS; level++) {
	assert(gkp_wa<T>(level).size() == fv1.size());
	assert(gkp_wa<T>(level).size() == fv2.size());
	assert(gkp_wb<T>(level).size() == gkp_x<T>(level).size()+1);
	T area2 = gkp_wb<T>(level).back() * f_center;
	// resabs = approximation to integral of abs(f)
	if (calcabsasc) resabs = fabs(area2);
	for (uint k=0; k<fv1.size(); k++) {
	  area2 += gkp_wa<T>(level)[k] * (fv1[k]+fv2[k]);
	  if (calcabsasc) resabs += gkp_wa<T>(level)[k] * (fabs(fv1[k]) + fabs(fv2[k]));
	}
	for (uint k=0; k<gkp_x<T>(level).size(); k++) {
	  const T abscissa = half_length * gkp_x<T>(level)[k];
	  const T fval1 = func(center - abscissa);
	  const T fval2 = func(center + abscissa);
	  const T fval = fval1 + fval2;
	  area2 += gkp_wb<T>(level)[k] * fval;
	  if (calcabsasc) resabs += gkp_wb<T>(level)[k] * (fabs(fval1) + fabs(fval2));
	  fv1.push_back(fval1);
	  fv2.push_back(fval2);
	  if (retfxmap) {
	    (*retfxmap)[center-abscissa] = fval1;
	    (*retfxmap)[center+abscissa] = fval2;
	  }
	}
#ifdef COUNTFEVAL
	nfeval+=gkp_x<T>(level).size()*2;
#endif
	if (calcabsasc) {
	  const T mean = area1*T(0.5);
	  // resasc = approximation to the integral of abs(f-mean) 
	  resasc = gkp_wb<T>(level).back() * fabs(f_center-mean);
	  for (uint k=0; k<gkp_wa<T>(level).size(); k++) {
	    resasc += gkp_wa<T>(level)[k] * (fabs(fv1[k]-mean) + fabs(fv2[k]-mean));
	  }
	  for (uint k=0; k<gkp_x<T>(level).size(); k++) {
	    resasc += gkp_wb<T>(level)[k] * (fabs(fv1[k]-mean) + fabs(fv2[k]-mean));
	  }
	  resasc *= abs_half_length ;
	  resabs *= abs_half_length;
	}
	area2 *= half_length;
	err = rescale_error (fabs(area2-area1), resabs, resasc) ;
	if (err < resasc) calcabsasc = false;

	integ_dbg2<<"at level "<<level<<" area2 = "<<area2<<" +- "<<err<<std::endl;

	//   test for convergence.
	if (err < epsabs || err < epsrel * fabs (area2)) {
	  reg.SetArea(area2,err);
	  return true;
	}
	area1 = area2;
      }

      // failed to converge 
      reg.SetArea(area1,err);

      integ_dbg2<<"failed to reach tolerance with highest-order GKP rule\n";
      return false;
    }

  template <class T>
    T FindPeriod(const std::vector<T>& v) 
    {
      uint minn = v.size()/2;
      uint n=1;
      while (n<=minn) n*=2;

      std::vector<T> dv(v.begin() + v.size()-n,v.end());
      assert(dv.size() == n);
      // Welch window:
      uint no2 = n/2;
      for(uint i=0;i<no2;i++) {
	T w = (no2-i)/no2;
	w = T(1)-w*w;
	dv[i] *= w;
	if (i>0) dv[n-i] *= w;
      }
      std::vector<std::complex<T> > fft;
      FFT(dv,fft);
      T maxpower = 0.;
      uint imax=0;
      for(uint i=1;i<fft.size();i++) {
	// start at i=1 to ignore DC power
	T power = norm(fft[i]);
	if (power > maxpower) {
	  imax = i; maxpower = power;
	}
      }
      assert(imax>0);
      if (imax == 1) return 0.;

      T dimax=0., sumpower=0.;
      assert(imax >= 2);
      for(uint i=imax-2;i<=imax+2;i++) {
	if (i >= fft.size()) continue;
	T power = norm(fft[i]);
	dimax += power*i;
	sumpower += power;
      }
      assert(sumpower > 0.);
      dimax /= sumpower;
      assert(dimax > 0.);
      return T(n)/dimax;
    }

  // Slightly different if T is complex:
  template <class T>
    T FindPeriod(const std::vector<std::complex<T> >& v) 
    {
      uint minn = v.size()/2;
      uint n=1;
      while (n<=minn) n*=2;

      std::vector<std::complex<T> > fft(n);
      for(uint i=0;i<n;i++) {
	const std::complex<T>& c = v[i+v.size()-n];
	fft[i] = std::complex<T>(real(c),imag(c));
      }
      // Welch window:
      uint no2 = n/2;
      for(uint i=0;i<no2;i++) {
	T w = (no2-i)/no2;
	w = T(1)-w*w;
	fft[i] *= w;
	if(i>0) fft[n-i] *= w;
      }
      FFT(fft);
      T maxpower = 0.;
      uint imax=0;
      std::vector<T> power(n/2+2,0.);
      for(uint i=1;i<=n/2;i++) {
	// start at i=1 to ignore DC power
	power[i] = norm(fft[i])+norm(fft[n-i]);
	if (power[i] > maxpower) {
	  imax = i; maxpower = power[i];
	}
      }
      if (imax == 1) return 0.;

      T dimax=0., sumpower=0.;
      for(uint i=imax-2;i<=imax+2;i++) {
	dimax += power[i]*i;
	sumpower += power[i];
      }
      dimax /= sumpower;
      return T(n)/dimax;
    }

  template <class UF>
    bool IntOsc (const UF& func, IntRegion<typename UF::result_type>& reg,
	const typename UF::result_type& period,
	const typename UF::result_type epsabs,
	const typename UF::result_type epsrel, 
	const typename UF::result_type extraprec=100.)
    // An integration algorithm appropriate for integrals of oscillating
    // functions from a to infinity.  
    // It integrates the first 20 periods worth with GKPOSC
    // Then it successively finds the integral for one wavelength
    // at a time and extrapolates the sequence of the sum of these 
    // integrals to get the overall sum.
    // At least 10 wavelengths are calculated before trying to extrapolate
    //
    // Returns: true if succeed, false if fail
    // The area and estimated error are returned as reg.Area() and reg.Err()
    //
    {
      typedef typename UF::result_type T;
      assert(reg.Alg() == OSCINFTWOPI || reg.Alg() == OSCINFONE);

      assert(epsabs >= 0.);
      assert(epsrel > 0.);
      if (extraprec > 2.e10) {
	std::cerr<<"Extra precision > 10^10 asked for in IntOsc.  Give up.\n";
	return false;
      }

      T x0 = reg.Left();
      T x1 = x0 + 5.*period;
      IntRegion<T> startreg(x0,x1,GKPOSC);
      if (!IntGKP(func,startreg,epsabs/extraprec,epsrel/extraprec)) {
	std::cerr<<"Integration of startreg failed in IntOsc\n";
	return false;
      }
      integ_dbg2<<"area of startreg = "<<startreg.Area()<<" +- "<<startreg.Err()<<std::endl;

      bool converged = false;
      const int MAXITER=64;
      SeriesLevin<T> areaseries;
      areaseries.AddTerm(startreg.Area(),startreg.Err());
      T interr = startreg.Err();
      T extraparea = startreg.Area();
      T extraperr = startreg.Err();
      for(int iter=1; iter<=MAXITER; iter++) {
	x0 = x1;
	x1 += period;
	IntRegion<T> oneperiod(x0,x1,GKPOSC);
	if (!IntGKP(func,oneperiod,epsabs/extraprec,epsrel/extraprec)) {
	  std::cerr<<"Integration of oneperiod failed in IntOsc\n";
	  return false;
	}
	integ_dbg2<<"area of period "<<iter<<" = "<<oneperiod.Area()<<" +- "<<oneperiod.Err()<<std::endl;
	areaseries.AddTerm(oneperiod.Area(),oneperiod.Err());
	interr += oneperiod.Err();
	areaseries.Extrapolate(&extraparea,&extraperr,0);
	extraperr *= 100.;  // It can be optimistic at times, up to 2 digits
	if (iter >= 5) {
	  T toler = std::max(epsabs,epsrel*fabs(extraparea));
	  integ_dbg2<<"area after "<<iter<<" oscillations = "<<extraparea;
	  integ_dbg2<<" +- "<<extraperr<<" (tol = "<<toler<<")\n";
	  if (extraperr < toler) {
	    converged = true;
	    break;
	  }
	  if (interr > toler) {
	    integ_dbg2<<"interr = "<<interr<<", epsabs = "<<epsabs<<", epsrel*area = "<<epsrel*fabs(extraparea)<<std::endl;
	    integ_dbg2<<"Integration errors already too big to ever converge\n";
	    integ_dbg2<<"Redo with extraprec = "<<extraprec*10<<std::endl;
	    return IntOsc(func,reg,period,epsabs,epsrel,extraprec*10);
	  }
	  if (areaseries.PastBest()) {
	    integ_dbg2<<"Series extrapolation is not going to improve.\n";
	    integ_dbg2<<"Redo with extraprec = "<<extraprec*10<<std::endl;
	    return IntOsc(func,reg,period,epsabs,epsrel,extraprec*10);
	  }
	} 
	else {
	  integ_dbg2<<"area after "<<iter<<" oscillations = "<<extraparea;
	  integ_dbg2<<" +- "<<extraperr<<std::endl;
	}
      }
      reg.SetArea(extraparea,extraperr);
      if (!converged) {
	integ_dbg2<<"Series extrapolation not good enough after MAXITER iterations\n";
	// Probably the period is not the best period to use
	// Find a better (longer) period for another pass:
	T newperiod = FindPeriod(areaseries.GetTerms());
	integ_dbg2<<"newperiod = "<<newperiod<<std::endl;
	if (newperiod == 0.0) {
	  integ_dbg2<<"There isn't a better period for this.  This one is just a kicker.\n";
	  newperiod = 100.;
	}
	integ_dbg2<<"Redo with period = "<<newperiod<<"*period (= "<<newperiod*period<<")\n";
	return IntOsc(func,reg,newperiod*period,epsabs,epsrel,extraprec);
      }
      return true;
    }

  template <class UF>
    bool IntGKP (const UF& func, IntRegion<typename UF::result_type>& reg,
	const typename UF::result_type epsabs,
	const typename UF::result_type epsrel,
	std::map<typename UF::result_type,typename UF::result_type>* retfxmap=0)
    // An adaptive integration algorithm which computes the integral of f
    // over the region reg.
    // First the non-adaptive GKP algorithm is tried.
    // If that is not accurate enough (according to the absolute and
    // relative accuracies, epsabs and epsrel),
    // the region is split in half, and each new region is integrated.
    // The routine continues by successively splitting the subregion
    // which gave the largest absolute error until the integral converges.
    //
    // Returns true is succeed, false if fail
    // The area and estimated error are returned as reg.Area() and reg.Err()
    // If desired, *retx and *retf return std::vectors of x,f(x) respectively
    // They only include the evaluations in the non-adaptive pass, so they
    // do not give an accurate estimate of the number of function evaluations.
    {
      typedef typename UF::result_type T;
      bool oscsplit = (reg.Alg() == GKPOSC);
      integ_dbg2<<"Start IntGKP with oscsplit = "<<oscsplit<<std::endl;
      assert(oscsplit || reg.Alg() == GKP);

      assert(epsabs >= 0.);
      assert(epsrel > 0.);

      // perform the first integration 
      bool done;
      std::map<T,T> fxmap;
      if (oscsplit) {
	done = IntGKPNA(func, reg, epsabs, epsrel, &fxmap);
	if (!done) reg.FindOscSplits(fxmap);
      } else {
	done = IntGKPNA(func, reg, epsabs, epsrel, &fxmap);
      }
      if (retfxmap) *retfxmap = fxmap;
      if (done) return true;

      integ_dbg2<<"In adaptive GKP, failed first pass... subdividing\n";
      integ_dbg2<<"Intial range = "<<reg.Left()<<".."<<reg.Right()<<std::endl;

      int roundoff_type1 = 0, error_type = 0;
      T roundoff_type2 = 0.;
      uint iteration = 1;

      std::priority_queue<IntRegion<T>,std::vector<IntRegion<T> > > allregions;
      allregions.push(reg);
      T finalarea = reg.Area();
      T finalerr = reg.Err();
      T tolerance= std::max(epsabs, epsrel * fabs(finalarea));
      assert(finalerr > tolerance);

      while(!error_type && finalerr > tolerance) {
	// Bisect the subinterval with the largest error estimate 
	integ_dbg2<<"Current answer = "<<finalarea<<" +- "<<finalerr<<"  (tol = "<<tolerance<<")\n";
	IntRegion<T> parent = allregions.top(); 
	allregions.pop();
	integ_dbg2<<"Subdividing largest error region "<<parent.Left()<<".."<<parent.Right()<<std::endl;
	integ_dbg2<<"parent area = "<<parent.Area()<<" +- "<<parent.Err()<<std::endl;
	std::vector<IntRegion<T> > children;
	parent.SubDivide(&children);
	// For "GKP", there are only two, but for GKPOSC, there is one 
	// for each oscillation in region

	// Try to do at least 3x better with the children
	T factor = 3*children.size()*finalerr/tolerance;
	T newepsabs = fabs(parent.Err()/factor);
	T newepsrel = newepsabs/fabs(parent.Area());
	integ_dbg2<<"New epsabs,rel = "<<newepsabs<<','<<newepsrel<<"  ("<<children.size()<<" children)\n";

	T newarea = T(0.0);
	T newerror = 0.0;
	for(uint i=0;i<children.size();i++) {
	  IntRegion<T>& child = children[i];
	  integ_dbg2<<"Integrating child "<<child.Left()<<".."<<child.Right()<<std::endl;
	  bool converged;
	  if (oscsplit) {
	    converged = IntGKPNA(func, child, newepsabs, newepsrel, &fxmap);
	    if (!converged) child.FindOscSplits(fxmap);
	  } else {
	    converged = IntGKPNA(func, child, newepsabs, newepsrel);
	  }
	  integ_dbg2<<"child ("<<i+1<<'/'<<children.size()<<") ";
	  if (converged) {integ_dbg2<<" converged.";}
	  else {integ_dbg2<<" failed.";}
	  integ_dbg2<<"  Area = "<<child.Area()<<" +- "<<child.Err()<<std::endl;

	  newarea += child.Area();
	  newerror += child.Err();
	}
	integ_dbg2<<"Compare: newerr = "<<newerror<<" to parent err = "<<parent.Err()<<std::endl;

	finalerr += (newerror - parent.Err());
	finalarea += newarea - parent.Area();

	T delta = parent.Area() - newarea;
	if (newerror <= parent.Err() && fabs (delta) <=  parent.Err() && newerror >= 0.99 * parent.Err()) {
	  integ_dbg2<<"roundoff type 1: delta/newarea = "<<fabs(delta)/fabs(newarea)<<", newerror/error = "<<newerror/parent.Err()<<std::endl;
	  roundoff_type1++;
	}
	if (iteration >= 10 && newerror > parent.Err() && fabs(delta) <= newerror-parent.Err()) {
	  integ_dbg2<<"roundoff type 2: newerror/error = "<<newerror/parent.Err()<<std::endl;
	  roundoff_type2+=std::min(newerror/parent.Err()-1.,T(1.));
	}

	tolerance = std::max(epsabs, epsrel * fabs(finalarea));
	if (finalerr > tolerance) {
	  if (roundoff_type1 >= 200) {
	    error_type = 1;	// round off error 
	    integ_dbg2<<"GKP: Round off error 1\n";
	  }
	  if (roundoff_type2 >= 200.) {
	    error_type = 2;	// round off error 
	    integ_dbg2<<"GKP: Round off error 2\n";
	  }
	  if (fabs((parent.Right()-parent.Left())/(reg.Right()-reg.Left())) < Epsilon<double>()) {
	    error_type = 3; // found singularity
	    integ_dbg2<<"GKP: Probable singularity\n";
	  }
	}
	for(uint i=0;i<children.size();i++) allregions.push(children[i]);
	iteration++;
      } 

      // Recalculate finalarea in case there are any slight rounding errors
      finalarea=0.; finalerr=0.;
      while (!allregions.empty()) {
	const IntRegion<T>& r=allregions.top();
	finalarea += r.Area();
	finalerr += r.Err();
	allregions.pop();
      }
      reg.SetArea(finalarea,finalerr);

      if (error_type == 1) {
	std::cerr << "Type 1 roundoff's = "<<roundoff_type1<<", Type 2 = "<<roundoff_type2<<std::endl;
	std::cerr << "Roundoff error 1 prevents tolerance from being achieved in IntGKP\n";
	return false;
      }
      else if (error_type == 2) {
	std::cerr << "Type 1 roundoff's = "<<roundoff_type1<<", Type 2 = "<<roundoff_type2<<std::endl;
	std::cerr << "Roundoff error 2 prevents tolerance from being achieved in IntGKP\n";
	return false;
      }
      else if (error_type == 3) {
	std::cerr << "Bad integrand behavior found in the integration interval in IntGKP\n";
	return false; 
      }
      else {
	return true;
      }
    }


  const double DEFABSERR = 1.e-15;
  const double DEFRELERR = 1.e-6;

  template <class UF>
    struct AuxFunc1   // f(1/x-1) for int(a..infinity)
    : public std::unary_function<typename UF::argument_type,
  typename UF::result_type> 
  {
    public:
      AuxFunc1(const UF& _f) : f(_f) {}
      typename UF::result_type 
	operator()(typename UF::argument_type x) const 
	{ return f(1./x-1.)/(x*x); }
    private:
      const UF& f;
  };

  template <class UF>
    struct AuxFunc1<UF> Aux1(UF uf) { return AuxFunc1<UF>(uf); }

  template <class UF>
    struct AuxFunc2  // f(1/x+1) for int(-infinity..b)
    : public std::unary_function<typename UF::argument_type,
  typename UF::result_type> 
  {
    public:
      AuxFunc2(const UF& _f) : f(_f) {}
      typename UF::result_type 
	operator()(typename UF::argument_type x) const 
	{ return f(1./x+1.)/(x*x); }
    private:
      const UF& f;
  };

  template <class UF>
    struct AuxFunc2<UF> Aux2(UF uf) { return AuxFunc2<UF>(uf); }

  template <class UF>
    struct AuxFunc3  // f(-x) for int(-infinity..b) for OSCINF
    : public std::unary_function<typename UF::argument_type,
  typename UF::result_type> 
  {
    public:
      AuxFunc3(const UF& _f) : f(_f) {}
      typename UF::result_type 
	operator()(typename UF::argument_type x) const 
	{ return f(-x); }
    private:
      const UF& f;
  };

  template <class UF>
    struct AuxFunc3<UF> Aux3(UF uf) { return AuxFunc3<UF>(uf); }

  template <class UF>
    typename UF::result_type int1d(const UF& func, 
	IntRegion<typename UF::result_type>& reg,
	const typename UF::result_type& abserr=DEFABSERR,
	const typename UF::result_type& relerr=DEFRELERR)
    {
      typedef typename UF::result_type T;

      integ_dbg2<<"start int1d: "<<reg.Left()<<".."<<reg.Right()<<std::endl;

      const T MOCKINF = 1.e10;
      const double PI = 3.141592653589793;
      const double TWOPI = 2.*PI;

      if ((reg.Left() <= -MOCKINF && reg.Right() > 0) ||
	  (reg.Right() >= MOCKINF && reg.Left() < 0)) { 
	reg.AddSplit(0);
      }

      if (reg.NSplit() > 0) {
	std::vector<IntRegion<T> > children;
	reg.SubDivide(&children);
	integ_dbg2<<"Subdivided into "<<children.size()<<" children\n";
	T answer=T();
	T err; 
	for(uint i=0;i<children.size();i++) {
	  IntRegion<T>& child = children[i];
	  integ_dbg2<<"i = "<<i<<": bounds = "<<child.Left()<<','<<child.Right()<<std::endl;
	  answer += int1d(func,child,abserr,relerr);
	  err += child.Err();
	  integ_dbg2<<"subint = "<<child.Area()<<" +- "<<child.Err()<<std::endl;
	}
	reg.SetArea(answer,err);
	return answer;
      }
      else {
	bool fail=false;
	T period=0.;
	switch(reg.Alg()) {
	  // Any additional OSCINF periods can go here and fall through
	  case OSCINFTWOPI : 
	    {
	      if(reg.Alg() == OSCINFTWOPI) period = TWOPI;
	    } // fall through...
	  case OSCINFONE : 
	    {
	      if(reg.Alg() == OSCINFONE) period = 1.;
	      integ_dbg2<<"oscinf, period = "<<period<<std::endl;
	      if (reg.Left() <= -MOCKINF) {
		integ_dbg2<<"left = -infinity, right = "<<reg.Right()<<std::endl;
		assert(reg.Right() <= 0.);
		IntRegion<T> modreg(-reg.Right(),-reg.Left(),reg.Alg(),reg.dbgout);
		if(!IntOsc(Aux3(func),modreg,period,abserr,relerr))
		{ fail = true; std::cerr<<"IntOsc failed\n"; }
		reg.SetArea(modreg.Area(),modreg.Err());
	      } else if (reg.Right() >= MOCKINF) {
		integ_dbg2<<"left = "<<reg.Left()<<", right = infinity\n";
		assert(reg.Left() >= 0.);
		if(!IntOsc(func,reg,period,abserr,relerr))
		{ fail = true; std::cerr<<"IntOsc failed\n"; }
	      } else {
		integ_dbg2<<"left = "<<reg.Left()<<", right = "<<reg.Right()<<std::endl;
		IntRegion<T> modreg(reg.Left(),reg.Right(),GKPOSC,reg.dbgout);
		if(!IntGKP(func,modreg,abserr,relerr))
		{ fail = true; std::cerr<<"IntGKP failed\n"; }
		reg.SetArea(modreg.Area(),modreg.Err());
	      }
	    } break;
	  case GKP : case GKPOSC : 
	    {
	      if (reg.Left() <= -MOCKINF) {
		integ_dbg2<<"left = -infinity, right = "<<reg.Right()<<std::endl;
		assert(reg.Right() <= 0.);
		IntRegion<T> modreg(1./(reg.Right()-1.),0.,reg.Alg(),reg.dbgout);
		if(!IntGKP(Aux2<UF>(func),modreg,abserr,relerr))
		{ fail = true; std::cerr<<"IntGKP failed\n"; }
		reg.SetArea(modreg.Area(),modreg.Err());
	      } else if (reg.Right() >= MOCKINF) {
		integ_dbg2<<"left = "<<reg.Left()<<", right = infinity\n";
		assert(reg.Left() >= 0.);
		IntRegion<T> modreg(0.,1./(reg.Left()+1.),reg.Alg(),reg.dbgout);
		if(!IntGKP(Aux1<UF>(func),modreg,abserr,relerr))
		{ fail = true; std::cerr<<"IntGKP failed\n"; }
		reg.SetArea(modreg.Area(),modreg.Err());
	      } else {
		integ_dbg2<<"left = "<<reg.Left()<<", right = "<<reg.Right()<<std::endl;
		if(!IntGKP(func,reg,abserr,relerr))
		{ fail = true; std::cerr<<"IntGKP failed\n"; }
	      }
	    } break;
	  default : 
	    {
	      myerror("Invalid algorithm");
	    } break;
	}
	if (fail) myerror("Integration failed"); 
	integ_dbg2<<"done int1d  answer = "<<reg.Area()<<" +- "<<reg.Err()<<std::endl;
	return reg.Area();
      }
    }

  template <class BF, class YREG>
    class Int2DAuxType : 
    public std::unary_function<typename BF::first_argument_type,typename BF::result_type> {
      public:
	Int2DAuxType(const BF& _func, const YREG& _yreg,
	    const typename BF::result_type& _abserr,
	    const typename BF::result_type& _relerr) :
	  func(_func),yreg(_yreg),abserr(_abserr),relerr(_relerr) {}

	typename BF::result_type operator()(typename BF::first_argument_type x) const {
	  typename YREG::result_type tempreg = yreg(x);
	  typename BF::result_type result = int1d(bind21(func,x),tempreg,abserr,relerr);
	  integ_dbg3<<"Evaluated int2dAux at x = "<<x<<": f = "<<result<<" +- "<<tempreg.Err()<<std::endl;
	  return result;
	} 

      private:
	const BF& func;
	const YREG& yreg;
	typename BF::result_type abserr,relerr;
    };

  template <class BF, class YREG>
    typename BF::result_type int2d(const BF& func,
	IntRegion<typename BF::result_type>& reg,
	const YREG& yreg, const typename BF::result_type& abserr=DEFABSERR,
	const typename BF::result_type& relerr=DEFRELERR)
    {
      integ_dbg2<<"Starting int2d: range = "<<reg.Left()<<".."<<reg.Right()<<std::endl;
      Int2DAuxType<BF,YREG> faux(func,yreg,abserr*1.e-3,relerr*1.e-3);
      typename BF::result_type answer = int1d(faux,reg,abserr,relerr);
      integ_dbg2<<"done int2d  answer = "<<answer<<" +- "<<reg.Err()<<std::endl;
      return answer;
    }

  template <class TF, class YREG, class ZREG>
    class Int3DAuxType : 
    public std::unary_function<typename TF::firstof3_argument_type,typename TF::result_type> {
      public:
	Int3DAuxType(const TF& _func, const YREG& _yreg, const ZREG& _zreg, 
	    const typename TF::result_type& _abserr, const typename TF::result_type& _relerr) :
	  func(_func),yreg(_yreg),zreg(_zreg),abserr(_abserr),relerr(_relerr) {}

	typename TF::result_type operator()(typename TF::firstof3_argument_type x) const {
	  typename YREG::result_type tempreg = yreg(x);
	  typename TF::result_type result = 
	    int2d(bind31(func,x),tempreg,bind21(zreg,x),abserr,relerr);
	  integ_dbg3<<"Evaluated int3dAux at x = "<<x<<": f = "<<result<<" +- "<<tempreg.Err()<<std::endl;
	  return result;
	}

      private:
	const TF& func;
	const YREG& yreg;
	const ZREG& zreg;
	typename TF::result_type abserr,relerr;
    };

  template <class TF, class YREG, class ZREG>
    typename TF::result_type int3d(const TF& func,
	IntRegion<typename TF::result_type>& reg, const YREG& yreg, const ZREG& zreg,
	const typename TF::result_type& abserr=DEFABSERR,
	const typename TF::result_type& relerr=DEFRELERR)
    {
      integ_dbg2<<"Starting int3d: range = "<<reg.Left()<<".."<<reg.Right()<<std::endl;
      Int3DAuxType<TF,YREG,ZREG> faux(func,yreg,zreg,abserr*1.e-3,relerr*1.e-3);
      typename TF::result_type answer = int1d(faux,reg,abserr,relerr);
      integ_dbg2<<"done int3d  answer = "<<answer<<" +- "<<reg.Err()<<std::endl;
      return answer;
    }

  // Helpers for constant regions for int2d, int3d:

  template <class T>
    struct ConstantReg1 : std::unary_function<T, IntRegion<T> >
    {
      ConstantReg1(T a,T b,INTALG alg=DEFALG) : ir(a,b,alg) {}
      IntRegion<T> operator()(T x) const {return ir;}
      IntRegion<T> ir;
    };

  template <class T>
    struct ConstantReg2 : std::binary_function<T, T, IntRegion<T> >
    {
      ConstantReg2(T a,T b,INTALG alg=DEFALG) : ir(a,b,alg) {}
      IntRegion<T> operator()(T x, T y) const {return ir;}
      IntRegion<T> ir;
    };

}

#endif
